{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><h1> Latent linear models</h1></center>\n",
    "\n",
    "# 1. Factor analysis (FA)\n",
    "One problem with mixture models is that they only use a single latent variable to generate the observations. In particular, each observation can only come from one of $K$ prototypes. An alternative is to use a vector of real-valued latent variables $\\vec{z_i} \\in R^L$. The simplest prior to use is a Gaussian\n",
    "$$\n",
    "p(\\vec{z}_i)=\\mathcal{N}(\\vec{z}_i|\\vec{\\mu}_0,\\Sigma_0)\n",
    "$$\n",
    "If the observations are also continuous, so $\\vec{x_i} \\in R^D$, we may use a Gaussian for the conditional distribution \n",
    "$$\n",
    "p(\\vec{x}_i|\\vec{z}_i,\\theta)=\\mathcal{N}(\\vec{x}_i|\\mathbf{W}\\vec{z}_i+\\vec{\\mu},\\mathbf{\\Psi})\n",
    "$$\n",
    "where $\\mathbf{W}$ is a $ D \\times L$ matrix,known as the **factor loading matrix**. And $\\Psi$ is a $D \\times D$ covariance matrix, which is always diagonal, because the aim of the factor analysis model is to \"force\" the latent $\\vec{z}_i$ to explain the correlation in $\\vec{x}_i$.\n",
    "\n",
    "## 1.1 FA is a low rank parameterization of an MVN\n",
    "FA can be thought of as a way of specifying a joint density model on $\\vec{x}$ using a small number of parameters. From the defination of FA, we can compute the induced marginal distribution as follows\n",
    "\\begin{align}\n",
    "p(\\vec{x_i}|\\theta) &= \\int \\mathcal{N}(\\vec{x}_i|\\mathbf{W}\\vec{z}_i+\\mu,\\Psi) \\mathcal{N}(\\vec{z}_i|\\vec{\\mu}_0,\\Sigma_0) d\\vec{z}_i \\\\\n",
    "                    &=\\mathcal{N}(\\vec{x}_i|\\mathbf{W}\\vec{\\mu}_0+\\vec{\\mu},\\Psi+\\mathbf{W}\\Sigma_0\\mathbf{W}^T) \n",
    "\\end{align}\n",
    "We can set $\\vec{\\mu}_0=\\mathbf{0}$ and $\\mathbf{\\Sigma_0}=\\mathbf{I}$ because we can always using linear transformation to do this. This means that the latent variables are just $L$ idd standard normal random variables.Then we get\n",
    "$$\n",
    "cov(\\vec{x}|\\theta)=\\Psi+\\mathbf{W}\\mathbf{W}^T\n",
    "$$\n",
    "This is just a low rank approximates of the covariance of an MVN.\n",
    "\n",
    "The FA can be also viewed as assuming that the data lives on a low dimensional linear manifold, which means\n",
    "$$\n",
    "\\mathbf{W}\\vec{z}_i=\\sum_{l=1}^L z_{il}\\vec{w}_l\n",
    "$$\n",
    "The data can be viewed as located in the subspace expanded by the vectors $\\lbrace \\vec{w}_1,\\vec{w}_2,\\ldots,\\vec{w}_L \\rbrace$.\n",
    "\n",
    "In the FA, we think the data lives on a low dimensional **linear manifold**. In reality, most data is better modeled by some form of low dimensional **curved manifold**, which is called **mixture of factor analysers**. We can approximate a curved mannifold by a piecewise linear manifold.\n",
    "\\begin{align}\n",
    "p(\\vec{x}_i|\\vec{z}_i,q_i=k,\\theta) &=\\mathcal{N}(\\vec{x}_i|\\mathbf{W}_k\\vec{z}_i+\\vec{\\mu}_k,\\Psi) \\\\\n",
    "                    p(\\vec{z}_i|\\theta) &=\\mathcal{N}(\\vec{z}_i|\\mathbf{0},\\mathbf{I}) \\\\\n",
    "                    p(q_i|\\pi)         &=Cat(q_i|\\pi)\n",
    "\\end{align}\n",
    "\n",
    "## 1.2 Inference of the lantent factors\n",
    "We can compute the posterior over the latent factors given the parameters using the Bayes rule.\n",
    "\\begin{align}\n",
    "p(\\vec{z}_i|\\vec{x}_i,\\theta) &\\propto p(\\vec{z}_i|\\theta)p(\\vec{x}_i|\\vec{z}_i,\\theta) \\\\\n",
    "                              &= \\mathcal{N}(\\vec{z}_i|\\vec{\\mu}_0,\\Sigma_0) \\mathcal{N}(\\vec{x}_i|\\mathbf{W}\\vec{z}_i+\\vec{\\mu},\\mathbf{\\Psi}) \\\\\n",
    "                              &=\\mathcal{N}(\\vec{z}_i|\\vec{m}_i,\\Sigma_i) \\\\\n",
    "                     \\Sigma_i &=\\left(\\Sigma_0^{-1}+\\mathbf{W}^T\\Psi^{-1}\\mathbf{W}\\right)^{-1} \\\\\n",
    "                    \\vec{m}_i &=\\Sigma_i \\left( \\mathbf{W}^T\\Psi^{-1}(\\vec{x}_i-\\vec{\\mu})+\\Sigma_0^{-1}\\vec{\\mu}_0 \\right )\n",
    "\\end{align}\n",
    "where $\\vec{m}_i$ are called the **latent scores**, or **latent factors**.\n",
    "\n",
    "## 1.3 Unidentifiability\n",
    "Just like with mixture models, FA is also unidentifiable. To see this ,suppose $\\mathbf{R}$ is an arbitrary orthogonal rotation matrix, which means $\\mathbf{RR}^T=\\mathbf{I}$. We can define $\\hat{\\mathbf{W}}=\\mathbf{WR}$;then the covariance of the marginal distribution $p(\\vec{x}|\\theta)$ is not changed. Therefore, we can't unique identify $\\mathbf{W}$\n",
    "\n",
    "To ensure a unique solution,we need to remove $\\frac{L(L-1)}{2}$ degrees of freedom,since that is the number of orthogonal rotation matrix of size $L \\times L$.There are some commonly used solutions\n",
    "\n",
    "* **Forcing $\\mathbf{W}$ to be orthonormal** The cleanest solution to the identifiability problem is to force $\\mathbf{W}$ to be orthonormal, which equals to that the columns of $\\mathbf{W}$ are orthonormal basis vector,  and to order the columns by decreasing variance of the corresponding latent factors.\n",
    "\\begin{align}\n",
    "\\mathbf{W}  &=\\lbrace \\vec{w}_1,\\ldots,\\vec{w}_L \\rbrace  \\\\\n",
    "\\vec{w}_i &\\perp \\vec{w}_j \\qquad  i,j=1,2,\\ldots,L\n",
    "\\end{align}\n",
    "* **Forcing $\\mathbf{W}$ to be lower triangular** One way to achieve identifiability, which is popular in the Bayesian community, is to ensure that the first visible feature is only generated by the first latent factor, the second visible feature is only generated by the first two latent factors.\n",
    "$$\n",
    "\\mathbf{W}= \\begin{bmatrix}\n",
    "w_{11} & 0 & 0 \\\\\n",
    "w_{21} & w_{22} & 0 \\\\\n",
    "w_{31} & w_{32} & w_{33} \\\\\n",
    "w_{41} & w_{42} & w_{43} \n",
    "\\end{bmatrix}  \n",
    "$$\n",
    "The first $L$ visible variables must be carefully chosen.\n",
    "* **Sparsity promoting priors on the weights** Instead of pre-specifying which entries in $\\mathbf{W}$ are zero, we can encourage the entries to be zero, using $\\mathcal{l}_1$ regularization, which is called sparse factor analysis.\n",
    "\n",
    "* **Use of non-Gaussian priors for the latent factors**. Replacing $p(\\vec{z}_i)$ with a non-Gaussian distribution can enable us to sometimes uniquely identify $\\mathbf{W}$ as well as the latent factors，which is known as ICA.\n",
    "\n",
    "# 2. Principal components analysis(PCA)\n",
    "Consider the FA model where we constrain $\\Psi=\\sigma^2\\mathbf{I}$ and $\\mathbf{W}$ to be orthonormal. As $ \\sigma^2 \\rightarrow  0$, this model reduces to classical (non-probabilistic) **principal components analysis(PCA)**. The version where $\\sigma^2 >0 $ is known as **probabilistic PCA(PPCA)**.\n",
    "\n",
    "## 2.1 Classical PCA\n",
    "The classical PCA can be viewed as a lossy compression method, which is similar to the  k-means method. But k-means method uses only a category scalar to represent a vector $\\vec{x}_i \\in R^D$, which losses lots of information. The classical PCA uses a vector of $\\vec{z}_i \\in R^L$ to represent a vector $\\vec{x}_i \\in R^D$. There are $L$ basis orthonormal vectors $\\vec{w}_l$, we need to minimize the average **reconstruction error**\n",
    "$$\n",
    "\\mathcal{J}(\\mathbf{W},\\mathbf{Z})=\\frac{1}{N}\\sum_{i=1}^N ||\\vec{x}_i-\\sum_{l=1}^L z_{il}\\vec{w}_l ||^2\n",
    "$$\n",
    "Equivalently,we can write this objective as follows\n",
    "$$\n",
    "\\mathcal{J}(\\mathbf{W},\\mathbf{Z})=||\\mathbf{X-WZ}^T||_F^2\n",
    "$$\n",
    "The optimal solution is obtained by setting $\\hat{\\mathbf{W}}=V_L$, where $V_L$ contains the $L$ eigenvectors with largest eigenvalues of the empirical covariance matrix $\\hat{\\Sigma}=\\frac{1}{N}\\sum_{i=1}^N (\\vec{x}_i-\\vec{\\mu})(\\vec{x}_i-\\vec{\\mu})^T$. Furthermore, the optimal low-dimensional encoding of the data is given by $\\vec{z}_i=\\mathbf{W}^T\\vec{x}_i$, which is an orthogonal projection of the data onto the column space spanned by the eigenvectors. This process is equal to finds the $l$ directions of maximal variance.This means that PCA can be “misled” by directions in which the variance is high merely because of the measurement scale. **It is therefore standard practice to standardize the data first.**\n",
    "\n",
    "## 2.2 Singular value decomposition(SVD)\n",
    "We have defined the solution of PCA in terms of eigenvectors of the covariance matrix. There is another way to obtain the solution ,based on the **singular value decomposition** of the matrix.\n",
    "\n",
    "In particular, any $N \\times D$ matrix $\\mathbf{X}$ can be decomposed as follows\n",
    "$$\n",
    "\\underbrace{\\mathbf{X}}_{N \\times D}=\\underbrace{\\mathbf{U}}_{N \\times N} \\quad \\underbrace{\\mathbf{S}}_{N \\times D} \\quad \\underbrace{\\mathbf{V}^T}_{D \\times D}\n",
    "$$\n",
    "where $U$ is an $N \\times N$ matrix which is orthornormal , $V$ is $D \\times D$ which is orthonormal, and $S$ is a $N \\times D$ matrix containing the $r=min(N,D)$ **singular values** $\\sigma_i \\geq 0$ on the main diagonal, with $0$ filling the rest of matrix.\n",
    "<img src=\"imgs/1.png\" alt=\"drawing\" width=\"500\"/>\n",
    "\n",
    "The connection between eigenvectors and singular vectors is the following. For and arbitrary real matrix $X$, we have\n",
    "\\begin{align}\n",
    "\\mathbf{X}^T\\mathbf{X} &=\\mathbf{VS}^T\\mathbf{U}^T\\mathbf{USV}^T \\\\\n",
    "                       &=\\mathbf{VS}^T\\mathbf{SV}^T\n",
    "\\end{align}\n",
    "which means the eigenvectors of $\\mathbf{X}^T\\mathbf{X}$ are equal to $\\mathbf{V}$, the right singular vectors of $\\mathbf{X}$, the eigenvalues of $\\mathbf{X}^T\\mathbf{X}$ are equal to the squared singular values of $\\mathbf{X}$\n",
    "\n",
    "From the SVD, we can represent a rank $r$ matrix as follows\n",
    "$$\n",
    "\\mathbf{X}=\\sigma_1\\vec{u}_1\\vec{v}_1^T+\\cdots+\\sigma_r\\vec{u}_r\\vec{v}_r^T\n",
    "$$\n",
    "If the singular values die off quickly, we can produce a rank $L$ approximation to the matrix as follows\n",
    "\\begin{align}\n",
    "\\mathbf{X} &=\\sigma_1\\vec{u}_1\\vec{v}_1^T+\\cdots+\\sigma_L\\vec{u}_L\\vec{v}_L^T  \\\\\n",
    "           &=\\mathbf{U}_{:,1:L}\\mathbf{S}_{1:L,1:L}\\mathbf{V}_{:,1:L}^T\n",
    "\\end{align}\n",
    "<img src=\"imgs/2.png\" alt=\"drawing\" width=\"500\"/>\n",
    "This is called **truncated SVD**, which is identity to the PCA.\n",
    "## 2.3 Probabilistic PCA\n",
    "Consider a factor analysis model in which $\\Psi=\\sigma^2\\mathbf{I}$ and $\\mathbf{W}$ is orthogonal. The observed  data log likelihood is given by (Tipping and Bishop 1999)\n",
    "\\begin{align}\n",
    "log\\,p(\\mathbf{X}|\\mathbf{W},\\sigma^2) &=-\\frac{N}{2}ln|\\mathbf{C}|-\\frac{1}{2}\\sum_{i=1}^{N}\\vec{x}_i^T\\mathbf{C}^{-1}\\vec{x}_i \\\\\n",
    "\\end{align}\n",
    "where $\\mathbf{C}=\\mathbf{W}\\mathbf{W}^T+\\sigma^2\\mathbf{I}$.\n",
    "The maxima of the log-likelihood are given by\n",
    "$$\n",
    "\\hat{\\mathbf{W}}=\\mathbf{V}(\\Lambda-\\sigma^2\\mathbf{I})^{\\frac{1}{2}}\\mathbf{R}\n",
    "$$\n",
    "where $\\mathbf{R}$ is an arbitrary $L \\times L$ orthogonal matrix, $\\mathbf{V}$ is the $D \\times L$ matrix whose columns are the first $L$ eigenvectors of the sample covariance matrix, and $\\Lambda$ is the corresponding diagonal matrix of eigenvalues.\n",
    "\n",
    "Although the usual way to fit a PCA model uses eigenvector methods, or the SVD, we can also use EM, which enable online learning of $\\mathbf{W}$. More details can be referred to this [paper](https://arxiv.org/pdf/1011.1745.pdf)\n",
    "## 2.4 PCA for categorical data\n",
    "When the data has the form $x_{ij} \\in \\lbrace 1,\\ldots, C \\rbrace,\\, i=1,\\ldots,N,\\, j=1,\\ldots,R$. We assume each $\\vec{x}_i$ is generated from a latent variable $\\vec{z}_i \\in R^L$, with a Gaussian prior, which is passed through the softmax function as follows\n",
    "\\begin{align}\n",
    "p(\\vec{z}_i) &=\\mathcal{N}(\\mathbf{0},\\mathbf{I}) \\\\\n",
    "p(\\vec{x}_i) &=\\prod_{j=1}^R Cat(x_{ij}|\\mathcal{S}\\left(\\mathbf{W}_j^T \\vec{z}_i+\\vec{w}_0\\right)\n",
    "\\end{align}\n",
    "where $\\mathbf{W}_j \\in R^{L \\times C} $ is the factor loading matrix for $x_j$ and $\\vec{w}_0$ is the offset term. We will call this **categorical PCA**.\n",
    "<img src=\"imgs/3.png\" alt=\"drawing\" width=\"500\"/>\n",
    "Left: 150 synthetic 16 dimensional bit vectors. Right: the 2d embedding learned by binary PCA, using variational EM. We have color coded points by the identity of the true “prototype” that generated them."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:tensorflow-cpu-1.11]",
   "language": "python",
   "name": "conda-env-tensorflow-cpu-1.11-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
